# Authors: Marco Portmann (marco.portmann@unifr.ch) with help for calculations by David Stadelmann 
# Please ask permission before using this skript in your projects
# Last change: 12.11.2015
#
#


#! Discrete Effects

#!2 Discrete Effect for each Variable holding the others constant

DiscreteEffects <- function(fit, evalspecial=NA, medianspecial=NA, baseline=F, printchange=F, skipFE = T){

require(alr3)

    # check variable name: Intercept and "=" not allowed
    IsProperVariable <- function(i){
        if( names(fit$coef)[i]!='Intercept' && regexpr("=", names(fit$coef)[i])==-1 & (i %in% fixedeffects & skipFE == T ) == F ){
        return(T)
        } else return(F)
      }


    # check whether variable with name vname has to be evaluated at
    # default value (=dval) or lower/upper exception value (=1, 2)

    evalexception <- function(vname, dval, lowup){
     if (length(evalspecial)>0  ){
      indx<- grep(vname, evalspecial[seq(from=1, to=length(evalspecial), by=3)])[1]
       if (length(indx)>0 & is.na(indx)==F){
           if (is.na( as.numeric(evalspecial[1+(indx-1)*3+lowup]))==F){
             return(as.numeric(evalspecial[1+(indx-1)*3+lowup]))
           }else{ stop("check evalspecial arguments")
                  return(dval)
                }
        }else{return(dval)}
     }else{return(dval)}
    }



    # construct Logit formula
    LogitText <- function(variables){
     txt <- ""
     for (i in 1:(length(variables)-1)){
       txt <- paste(txt, names(variables)[i], "*", variables[i], "+", sep=" ")
      }
     txt <- paste(txt, names(variables)[length(variables)], "*", variables[length(variables)],  sep=" ")

     txt <- paste( "exp(", txt, ")/ (1+ exp(",txt, "))", sep="")
     return(txt)
    }



    # compute discrete change and variance for variable i

    ComputeEffect <- function(i){

    # computes discrete effect, first variable changes, second is fixed
    
        # 1. take 1st and 3rd quartile

        firstQ <- summary(fit$x[, i - ShiftIntercept])[2]
        thirdQ <- summary(fit$x[, i - ShiftIntercept])[5]

        # 2. check whether it is a dummy and firstQ=thirdQ
        if(firstQ == thirdQ){
          if(thirdQ==0){
            thirdQ <- thirdQ + 1}
          else{
            firstQ <- firstQ - 1}
        }

        # 3. check whether an exception is defined
        evalData1 <- evalData
        evalData1[i] <- evalexception(names(evalData)[i], firstQ, 1)

        

        evalData2 <- evalData
       
        evalData2[i] <- evalexception(names(evalData)[i], thirdQ, 2)
     
        # replace interaction terms
        if (length(items)>0){
            # j in [list of all interaction terms which include variable i]
            for (j in unique(items[items[,1]==i,2])  ){
                # list of all variables in interaction term j
                tmp <- items[items[,2]==j,1]
                # replace the interaction term by the product of evalData1[tmp], i.e. by the product of all variables in the interaction
                if (length(tmp)>0){
                    evalData1[j] <- prod(evalData1[tmp])
                    evalData2[j] <- prod(evalData2[tmp])
                }         
            }
        }  
        
        
        # replace interaction terms' names
        for (k in iterms){

            tmpl <- items[ items[,2]==k,1]
            tstr <- ''
            for (j in 1:length(tmpl)){
              tstr <- paste(tstr, tmpl[[j]], sep='x')
            }
            colnames(RobVar)[k] <- tstr
            rownames(RobVar)[k] <- tstr
            names(fit$coefficients)[k] <-tstr
            names(evalData1)[k] <- tstr
            names(evalData2)[k] <- tstr
    }


        texE <- paste(LogitText(evalData2), "-", LogitText(evalData1) ,sep="")       
       # print(texE)
        discbE1 <- deltaMethod(coef(fit), vcov.=RobVar, g=texE)
        rownames(discbE1)[1] <- paste("DE",names(fit$coef)[i],sep="")

        if (printchange){
           discbE1 <- cbind(discbE1, LowerValue = evalData1[i], UpperValue = evalData2[i], Difference =  evalData2[i] - evalData1[i])
        }
    return(discbE1)

    }


# <-----------------------------------------------------------------------------
# (function starts here)

    RobVar <- fit$var
    if (is.null(RobVar)) RobVar <- vcov(fit)

    ShiftIntercept <- 0
    if (   length(grep("^Intercept$",  names(fit$coef)))>0 | length(grep("^\\(Intercept\\)$", names(fit$coef)))>0 ) ShiftIntercept <- 1

#    if (ShiftIntercept == 0) warning("No intercept found.")


    # list of interaction terms
    iterms <- grep('\\*', names(fit$coef))

    # list all variables which appear in an interaction    
    items <- matrix(0, nrow=0, ncol=2)

    # list all variables which appear in an interaction
    for(i in 1:length(coef(fit))){
      itemsinteractions <- unique(c( grep(paste(names(fit$coef)[i], '\\*', sep=' '), names(fit$coef)), grep(paste('\\*', names(fit$coef)[i] , sep=' '), names(fit$coef)) ))
      if (length(itemsinteractions)>0){
        items <-  rbind(items,
                cbind(i, (itemsinteractions)) )}
    }


    #        items <<- items
    #        iterms <<- iterms
    #        varx <<- colnames(RobVar)


    # vector with medians (interaction-variables are corrected later)
    if (ShiftIntercept>0){
    evalData <- apply(cbind(Intercept=rep(1,dim(fit$x)[1]),
                            fit$x), 2, median)
    } else{
    evalData <- apply(fit$x, 2, median)    
    
    }

    # replace names with =  AND SET VALUES TO 0 !!!!!
    
    fixedeffects <- grep("=", names(fit$coef))
    colnames(RobVar) <- gsub("=", "", colnames(RobVar))
    rownames(RobVar) <- gsub("=", "", rownames(RobVar))
    names(fit$coefficients) <- gsub("=", "", names(fit$coefficients))
    names(evalData) <- gsub("=", "", names(evalData))

    evalData[fixedeffects] <- 0


    


    #change evalData as given in medianspecial
    if (length(medianspecial)>1){
     evalData[medianspecial[seq(from=1, to=length(medianspecial), by=2)]] <- as.numeric(medianspecial[seq(from=2, to=length(medianspecial), by=2)])
    }


    # correct medians of interaction terms
    for (i in 1:length(iterms)){
       tmp <- items[items[,2]==iterms[i],1]
       if (length(tmp)>0)
        evalData[iterms[i]] <- prod(evalData[tmp])
    }

    # compute DE for variables
    # loop: drop interaction terms, intercept and variables with "=" in their name
    for (i in 1:length(coef(fit))){
     if ( IsProperVariable(i) & (i %in% iterms)==F){
            if ('tableDE' %in% ls()){
            tableDE <- rbind(tableDE, ComputeEffect(i))
            }else{
            tableDE <- ComputeEffect(i)
            }
     }
    }



    if (baseline==T & ShiftIntercept>0){
    
      texE <- " exp( Intercept * 1) / ( 1 + exp( Intercept * 1)) "
      discbE1 <- deltaMethod(coef(fit), vcov.=RobVar, g=texE)  
      rownames(discbE1)[1] <- paste("DE",names(fit$coef)[1],sep="")
      
     if (printchange){
       discbE1 <- cbind(discbE1, LowerValue = NA, UpperValue = NA, Difference =  NA)
     }
      
      tableDE <- rbind(discbE1, tableDE)
    }

    pvalue <- 2*(1-pnorm(abs(data.matrix(tableDE['Estimate']/tableDE['SE']))))
    colnames(pvalue)[1] <- 'pvalue'
    tableDE <- cbind(tableDE, pvalue)
class(tableDE) <- c(class(tableDE), "DiscreteEffects")

return(tableDE)
}


################################################################################
################################################################################

#!2 Changing several variables at once


# evalspecial  = c("Variable 1", lower value, upper value, "Variable 2", lower value, upper value, ....)
#                c("Variable 1", NA, NA, "Variable 2", lower value, upper value, ....)
#                  >>> NA sets value to 1th/3th quantile
#


DiscreteEffect <- function(fit, evalspecial, medianspecial=NA){

require(alr3)

    # check variable name: Intercept and "=" not allowed
    IsProperVariable <- function(i){
        if( names(fit$coef)[i]!='Intercept' && regexpr("=", names(fit$coef)[i])==-1 ){
        return(T)
        } else return(F)
      }


    # check whether variable with name vname has to be evaluated at
    # default value (=dval) or lower/upper exception value (=1, 2)

    evalexception <- function(vname, dval, lowup){
     if (length(evalspecial)>0  ){
      indx<- grep(vname, evalspecial[seq(from=1, to=length(evalspecial), by=3)])[1]
       if (length(indx)>0 & is.na(indx)==F){
           if (is.na( as.numeric(evalspecial[1+(indx-1)*3+lowup]))==F){
             return(as.numeric(evalspecial[1+(indx-1)*3+lowup]))
           }else{ stop("check evalspecial arguments")
                  return(dval)
                }
        }else{return(dval)}
     }else{return(dval)}
    }





    # construct Logit formula
    LogitText <- function(variables){
     txt <- ""
     for (i in 1:(length(variables)-1)){
       txt <- paste(txt, names(variables)[i], "*", variables[i], "+", sep=" ")
      }
     txt <- paste(txt, names(variables)[length(variables)], "*", variables[length(variables)],  sep=" ")

     txt <- paste( "exp(", txt, ")/ (1+ exp(",txt, "))", sep="")
     return(txt)
    }

    SetChangeVariables <- function(evalData, evalNames, evalValues, Quartile){

       evalValues <- as.numeric(evalValues) 

      # values given by evalValues
      if (length(evalNames[is.na(evalValues)==F])>0)
        evalData[evalNames[is.na(evalValues)==F]] <- evalValues[is.na(evalValues)==F]

      # if evalValues is NA use a quartile (for one value)
      if (length(evalNames[is.na(evalValues)])==1)
        evalData[evalNames[is.na(evalValues)]] <- quantile(fit$x[,evalNames[is.na(evalValues)]  ], 0.25*Quartile)
    
      # if evalValues is NA use a quartile (for vectors)
      if (length(evalNames[is.na(evalValues)])>1)
        evalData[evalNames[is.na(evalValues)]] <- apply(fit$x[,evalNames[is.na(evalValues)]  ], 2, quantile, 0.25*Quartile)
    

      # correct medians of interaction terms
      for (i in 1:length(iterms)){
        tmp <- items[items[,2]==iterms[i],1]
        if (length(tmp)>0)
          evalData[iterms[i]] <- prod(evalData[tmp])
      }
   
    return(evalData) 
    }

# <-----------------------------------------------------------------------------
# (function starts here)

    RobVar <- fit$var
    if (is.null(RobVar)) RobVar <- vcov(fit)

    ShiftIntercept <- 0
    if (   length(grep("^Intercept$",  names(fit$coef)))>0 | length(grep("^\\(Intercept\\)$", names(fit$coef)))>0 ) ShiftIntercept <- 1

#    if (ShiftIntercept == 0) warning("No intercept found.")



    # list of interaction terms
    iterms <- grep('\\*', names(fit$coef))
    # list all variables which appear in an interaction    
    items <- matrix(0, nrow=0, ncol=2)

    # list all variables which appear in an interaction
    for(i in 1:length(coef(fit))){
      itemsinteractions <- unique(c( grep(paste(names(fit$coef)[i], '\\*', sep=' '), names(fit$coef)), grep(paste('\\*', names(fit$coef)[i] , sep=' '), names(fit$coef)) ))
      if (length(itemsinteractions)>0){
        items <-  rbind(items,
                cbind(i, (itemsinteractions)) )}
    }

  # 1. Use median and medianspecial data to set values of evalData.
  #    Copy evalData to evalData1and evalData2. Overwrite all changed
  #    variables. 

    # vector with medians (interaction-variables are corrected later)
    if (ShiftIntercept>0){
    evalData <- apply(cbind(Intercept=rep(1,dim(fit$x)[1]),
                            fit$x), 2, median)
    } else{
    evalData <- apply(fit$x, 2, median)    
    
    }


    # replace names with =  AND SET VALUES TO 0 !!!!!
    
    fixedeffects <- grep("=", names(fit$coef))
    colnames(RobVar) <- gsub("=", "", colnames(RobVar))
    rownames(RobVar) <- gsub("=", "", rownames(RobVar))
    names(fit$coefficients) <- gsub("=", "", names(fit$coefficients))
    names(evalData) <- gsub("=", "", names(evalData))

    evalData[fixedeffects] <- 0





    #change evalData as given in medianspecial
    if (length(medianspecial)>1){
     evalData[medianspecial[seq(from=1, to=length(medianspecial), by=2)]] <- as.numeric(medianspecial[seq(from=2, to=length(medianspecial), by=2)])
    }



  # 2. Set lower and upper values for change variables and fix interaction values.
    evalData1 <- SetChangeVariables(evalData, 
                                    evalspecial[seq(from=1, to=length(evalspecial), by=3)],
                                    evalspecial[seq(from=2, to=length(evalspecial), by=3)],
                                    1)

    evalData2 <- SetChangeVariables(evalData, 
                                    evalspecial[seq(from=1, to=length(evalspecial), by=3)],
                                    evalspecial[seq(from=3, to=length(evalspecial), by=3)],
                                    3)
     
  # replace interaction terms' names
        for (k in iterms){

            tmpl <- items[ items[,2]==k,1]
            tstr <- ''
            for (j in 1:length(tmpl)){
              tstr <- paste(tstr, tmpl[[j]], sep='x')
            }
            colnames(RobVar)[k] <- tstr
            rownames(RobVar)[k] <- tstr
            names(fit$coefficients)[k] <-tstr
            names(evalData1)[k] <- tstr
            names(evalData2)[k] <- tstr
    }

   
   
    texE <- paste(LogitText(evalData2), "-", LogitText(evalData1) ,sep="")
    tableDE <- deltaMethod(coef(fit), vcov.=RobVar, g=texE)
    rownames(tableDE) <- "Discrete Change"


    pvalue <- 2*(1-pnorm(abs(data.matrix(tableDE['Estimate']/tableDE['SE']))))
    colnames(pvalue)[1] <- 'pvalue'
    tableDE <- cbind(tableDE, pvalue)
class(tableDE) <- c(class(tableDE), "DiscreteEffects")
return(tableDE)
}






#!2 evalspecial, set interactions manually


DiscreteEffectIntA <- function(fit, evalspecial, medianspecial=NA){

require(alr3)

    # check variable name: Intercept and "=" not allowed
    IsProperVariable <- function(i){
        if( names(fit$coef)[i]!='Intercept' && regexpr("=", names(fit$coef)[i])==-1 ){
        return(T)
        } else return(F)
      }


    # check whether variable with name vname has to be evaluated at
    # default value (=dval) or lower/upper exception value (=1, 2)

    evalexception <- function(vname, dval, lowup){
     if (length(evalspecial)>0  ){
      indx<- grep(vname, evalspecial[seq(from=1, to=length(evalspecial), by=3)])[1]
       if (length(indx)>0 & is.na(indx)==F){
           if (is.na( as.numeric(evalspecial[1+(indx-1)*3+lowup]))==F){
             return(as.numeric(evalspecial[1+(indx-1)*3+lowup]))
           }else{ stop("check evalspecial arguments")
                  return(dval)
                }
        }else{return(dval)}
     }else{return(dval)}
    }





    # construct Logit formula
    LogitText <- function(variables){
     txt <- ""
     for (i in 1:(length(variables)-1)){
       txt <- paste(txt, names(variables)[i], "*", variables[i], "+", sep=" ")
      }
     txt <- paste(txt, names(variables)[length(variables)], "*", variables[length(variables)],  sep=" ")

     txt <- paste( "exp(", txt, ")/ (1+ exp(",txt, "))", sep="")
     return(txt)
    }

    SetChangeVariables <- function(evalData, evalNames, evalValues, Quartile){

       evalValues <- as.numeric(evalValues) 

      # values given by evalValues
      if (length(evalNames[is.na(evalValues)==F])>0)
        evalData[evalNames[is.na(evalValues)==F]] <- evalValues[is.na(evalValues)==F]

      # if evalValues is NA use a quartile (for one value)
      if (length(evalNames[is.na(evalValues)])==1)
        evalData[evalNames[is.na(evalValues)]] <- quantile(fit$x[,evalNames[is.na(evalValues)]  ], 0.25*Quartile)
    
      # if evalValues is NA use a quartile (for vectors)
      if (length(evalNames[is.na(evalValues)])>1)
        evalData[evalNames[is.na(evalValues)]] <- apply(fit$x[,evalNames[is.na(evalValues)]  ], 2, quantile, 0.25*Quartile)
    

      # correct medians of interaction terms
      for (i in 1:length(iterms)){
        tmp <- items[items[,2]==iterms[i],1]
        if (length(tmp)>0)
          evalData[iterms[i]] <- prod(evalData[tmp])
      }
   
    return(evalData) 
    }

# <-----------------------------------------------------------------------------
# (function starts here)

    RobVar <- fit$var
    if (is.null(RobVar)) RobVar <- vcov(fit)

    ShiftIntercept <- 0
    if (   length(grep("^Intercept$",  names(fit$coef)))>0 | length(grep("^\\(Intercept\\)$", names(fit$coef)))>0 ) ShiftIntercept <- 1

#    if (ShiftIntercept == 0) warning("No intercept found.")
    replist <- data.frame(orig = as.character(), repl = as.character())

    
    slist <- c(evalspecial[seq(from=1, to=length(evalspecial), by=3)], medianspecial[seq(from=1, to=length(medianspecial), by=2)])
    
    slist <- slist[grep("\\*", slist)]
    slist <- slist[is.na(slist) == F]
    if (length(slist) > 0){
      replist <- data.frame(orig = slist, repl = paste("RrRVarXx",1:length(slist), sep = "_" ) )

      for (i in 1:length(replist)){      
        sltext <- paste("^", gsub("\\*", "\\\\*", replist[i, "orig"]), "$" , sep = "" )
        evalspecial <- gsub(sltext, replist[i, "repl"], evalspecial)
        medianspecial <- gsub(sltext, replist[i, "repl"], medianspecial)       
        names(fit$coef) <- gsub(sltext, replist[i, "repl"], names(fit$coef))
        colnames(RobVar) <- gsub(sltext, replist[i, "repl"], colnames(RobVar))
        rownames(RobVar) <- gsub(sltext, replist[i, "repl"], rownames(RobVar))
        names(fit$coefficients) <- gsub(sltext, replist[i, "repl"], names(fit$coefficients))
      }
    }

    # list of interaction terms
    iterms <- grep('\\*', names(fit$coef))
    # list all variables which appear in an interaction    
    items <- matrix(0, nrow=0, ncol=2)

    # list all variables which appear in an interaction
    for(i in 1:length(coef(fit))){
      itemsinteractions <- unique(c( grep(paste(names(fit$coef)[i], '\\*', sep=' '), names(fit$coef)), grep(paste('\\*', names(fit$coef)[i] , sep=' '), names(fit$coef)) ))
      if (length(itemsinteractions)>0){
        items <-  rbind(items,
                cbind(i, (itemsinteractions)) )}
    }

  # 1. Use median and medianspecial data to set values of evalData.
  #    Copy evalData to evalData1and evalData2. Overwrite all changed
  #    variables. 

    # vector with medians (interaction-variables are corrected later)
    if (ShiftIntercept>0){
    evalData <- apply(cbind(Intercept=rep(1,dim(fit$x)[1]),
                            fit$x), 2, median)
    } else{
    evalData <- apply(fit$x, 2, median)    
    
    }



    #change evalData as given in medianspecial
    if (length(medianspecial)>1){
     evalData[medianspecial[seq(from=1, to=length(medianspecial), by=2)]] <- as.numeric(medianspecial[seq(from=2, to=length(medianspecial), by=2)])
    }



  # 2. Set lower and upper values for change variables and fix interaction values.
    evalData1 <- SetChangeVariables(evalData, 
                                    evalspecial[seq(from=1, to=length(evalspecial), by=3)],
                                    evalspecial[seq(from=2, to=length(evalspecial), by=3)],
                                    1)

    evalData2 <- SetChangeVariables(evalData, 
                                    evalspecial[seq(from=1, to=length(evalspecial), by=3)],
                                    evalspecial[seq(from=3, to=length(evalspecial), by=3)],
                                    3)
     
  # replace interaction terms' names
        for (k in iterms){

            tmpl <- items[ items[,2]==k,1]
            tstr <- ''
            for (j in 1:length(tmpl)){
              tstr <- paste(tstr, tmpl[[j]], sep='x')
            }
            colnames(RobVar)[k] <- tstr
            rownames(RobVar)[k] <- tstr
            names(fit$coefficients)[k] <-tstr
            names(evalData1)[k] <- tstr
            names(evalData2)[k] <- tstr
    }

   
   
    texE <- paste(LogitText(evalData2), "-", LogitText(evalData1) ,sep="")
    tableDE <- deltaMethod(coef(fit), vcov.=RobVar, g=texE)
    rownames(tableDE) <- "Discrete Change"


    pvalue <- 2*(1-pnorm(abs(data.matrix(tableDE['Estimate']/tableDE['SE']))))
    colnames(pvalue)[1] <- 'pvalue'
    tableDE <- cbind(tableDE, pvalue)
class(tableDE) <- c(class(tableDE), "DiscreteEffects")
return(tableDE)
}



#!2 Difference in Difference




DiscreteEffectDiD <- function(fit, group, treatment, medianspecial=NA){

require(alr3)

    # check variable name: Intercept and "=" not allowed
    IsProperVariable <- function(i){
        if( names(fit$coef)[i]!='Intercept' && regexpr("=", names(fit$coef)[i])==-1 ){
        return(T)
        } else return(F)
      }


    # check whether variable with name vname has to be evaluated at
    # default value (=dval) or lower/upper exception value (=1, 2)

    evalexception <- function(vname, dval, lowup){
     if (length(evalspecial)>0  ){
      indx<- grep(vname, evalspecial[seq(from=1, to=length(evalspecial), by=3)])[1]
       if (length(indx)>0 & is.na(indx)==F){
           if (is.na( as.numeric(evalspecial[1+(indx-1)*3+lowup]))==F){
             return(as.numeric(evalspecial[1+(indx-1)*3+lowup]))
           }else{ stop("check evalspecial arguments")
                  return(dval)
                }
        }else{return(dval)}
     }else{return(dval)}
    }





    # construct Logit formula
    LogitText <- function(variables){
     txt <- ""
     for (i in 1:(length(variables)-1)){
       txt <- paste(txt, names(variables)[i], "*", variables[i], "+", sep=" ")
      }
     txt <- paste(txt, names(variables)[length(variables)], "*", variables[length(variables)],  sep=" ")

     txt <- paste( "exp(", txt, ")/ (1+ exp(",txt, "))", sep="")
     return(txt)
    }

    SetChangeVariables <- function(evalData, evalNames, evalValues, Quartile){

       evalValues <- as.numeric(evalValues)

      # values given by evalValues
      if (length(evalNames[is.na(evalValues)==F])>0)
        evalData[evalNames[is.na(evalValues)==F]] <- evalValues[is.na(evalValues)==F]

      # if evalValues is NA use a quartile (for one value)
      if (length(evalNames[is.na(evalValues)])==1)
        evalData[evalNames[is.na(evalValues)]] <- quantile(fit$x[,evalNames[is.na(evalValues)]  ], 0.25*Quartile)

      if (length(evalNames[is.na(evalValues)])>1)
      stop("Discrete Effects DiD Funktion �berpr�fen: Fehlende Werte werden auf Quartile = 1 gesetzt")

      # if evalValues is NA use a quartile (for vectors)
      if (length(evalNames[is.na(evalValues)])>1)
        evalData[evalNames[is.na(evalValues)]] <- apply(fit$x[,evalNames[is.na(evalValues)]  ], 2, quantile, 0.25*Quartile)


      # correct medians of interaction terms
      for (i in 1:length(iterms)){
        tmp <- items[items[,2]==iterms[i],1]
        if (length(tmp)>0)
          evalData[iterms[i]] <- prod(evalData[tmp])
      }

    return(evalData)
    }

# <-----------------------------------------------------------------------------
# (function starts here)



    RobVar <- fit$var
    if (is.null(RobVar)) RobVar <- vcov(fit)

    ShiftIntercept <- 0
    if (   length(grep("^Intercept$",  names(fit$coef)))>0 | length(grep("^\\(Intercept\\)$", names(fit$coef)))>0 ) ShiftIntercept <- 1

#    if (ShiftIntercept == 0) warning("No intercept found.")



    # list of interaction terms
    iterms <- grep('\\*', names(fit$coef))
    # list all variables which appear in an interaction
    items <- matrix(0, nrow=0, ncol=2)

    # list all variables which appear in an interaction
    for(i in 1:length(coef(fit))){
      itemsinteractions <- unique(c( grep(paste(names(fit$coef)[i], '\\*', sep=' '), names(fit$coef)), grep(paste('\\*', names(fit$coef)[i] , sep=' '), names(fit$coef)) ))
      if (length(itemsinteractions)>0){
        items <-  rbind(items,
                cbind(i, (itemsinteractions)) )}
    }

  # 1. Use median and medianspecial data to set values of evalData.
  #    Copy evalData to evalData1and evalData2. Overwrite all changed
  #    variables.

    # vector with medians (interaction-variables are corrected later)
    if (ShiftIntercept>0){
    evalData <- apply(cbind(Intercept=rep(1,dim(fit$x)[1]),
                            fit$x), 2, median)
    } else{
    evalData <- apply(fit$x, 2, median)

    }


    # replace names with =  AND SET VALUES TO 0 !!!!!
    
    fixedeffects <- grep("=", names(fit$coef))
    colnames(RobVar) <- gsub("=", "", colnames(RobVar))
    rownames(RobVar) <- gsub("=", "", rownames(RobVar))
    names(fit$coefficients) <- gsub("=", "", names(fit$coefficients))
    names(evalData) <- gsub("=", "", names(evalData))

    evalData[fixedeffects] <- 0



    #change evalData as given in medianspecial
    if (length(medianspecial)>1){
     evalData[medianspecial[seq(from=1, to=length(medianspecial), by=2)]] <- as.numeric(medianspecial[seq(from=2, to=length(medianspecial), by=2)])
    }




    gbefore   <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(1, 0), 1)

    gafter   <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(1, 1), 1)

    cbefore   <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(0, 0), 1)

    cafter    <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(0, 1), 1)


  # replace interaction terms' names
        for (k in iterms){

            tmpl <- items[ items[,2]==k,1]
            tstr <- ''
            for (j in 1:length(tmpl)){
              tstr <- paste(tstr, tmpl[[j]], sep='x')
            }
            colnames(RobVar)[k] <- tstr
            rownames(RobVar)[k] <- tstr
            names(fit$coefficients)[k] <-tstr
            names(gbefore)[k] <- tstr
            names(gafter)[k] <- tstr
            names(cbefore)[k] <- tstr
            names(cafter)[k] <- tstr
    }



    texE <- paste( "(", LogitText(gafter), "-", LogitText(gbefore), ") - (", LogitText(cafter), "-", LogitText(cbefore), ")", sep="")
    tableDE <- deltaMethod(coef(fit), vcov.=RobVar, g=texE)
    rownames(tableDE) <- "Discrete Change"


    pvalue <- 2*(1-pnorm(abs(data.matrix(tableDE['Estimate']/tableDE['SE']))))
    colnames(pvalue)[1] <- 'pvalue'
    tableDE <- cbind(tableDE, pvalue)
class(tableDE) <- c(class(tableDE), "DiscreteEffects")
return(tableDE)
}





#!2 Discrete Effects Treatment







DiscreteEffectTreatment <- function(fit, group, treatment, medianspecial=NA){
#Implementiert f�r "Difference-in-Difference" nach
#Puhani, P. A. The treatment effect, the cross difference, and the interaction term in nonlinear 'difference-in-differences' models Economics Letters, 2012, 115, 85-87
#Standardfehler der "Difference-in-Difference" wird mit Delta-Methode berechnet

require(alr3)

    # check variable name: Intercept and "=" not allowed
    IsProperVariable <- function(i){
        if( names(fit$coef)[i]!='Intercept' && regexpr("=", names(fit$coef)[i])==-1 ){
        return(call)
        } else return(F)
      }


    # check whether variable with name vname has to be evaluated at
    # default value (=dval) or lower/upper exception value (=1, 2)

    evalexception <- function(vname, dval, lowup){
     if (length(evalspecial)>0  ){
      indx<- grep(vname, evalspecial[seq(from=1, to=length(evalspecial), by=3)])[1]
       if (length(indx)>0 & is.na(indx)==F){
           if (is.na( as.numeric(evalspecial[1+(indx-1)*3+lowup]))==F){
             return(as.numeric(evalspecial[1+(indx-1)*3+lowup]))
           }else{ stop("check evalspecial arguments")
                  return(dval)
                }
        }else{return(dval)}
     }else{return(dval)}
    }





    # construct Logit formula
    LogitText <- function(variables){
     txt <- ""
     for (i in 1:(length(variables)-1)){
       txt <- paste(txt, names(variables)[i], "*", variables[i], "+", sep=" ")
      }
     txt <- paste(txt, names(variables)[length(variables)], "*", variables[length(variables)],  sep=" ")

     txt <- paste( "exp(", txt, ")/ (1+ exp(",txt, "))", sep="")
     return(txt)
    }

    SetChangeVariables <- function(evalData, evalNames, evalValues, Quartile){

       evalValues <- as.numeric(evalValues)

      # values given by evalValues
      if (length(evalNames[is.na(evalValues)==F])>0)
        evalData[evalNames[is.na(evalValues)==F]] <- evalValues[is.na(evalValues)==F]

      # if evalValues is NA use a quartile (for one value)
      if (length(evalNames[is.na(evalValues)])==1)
        evalData[evalNames[is.na(evalValues)]] <- quantile(fit$x[,evalNames[is.na(evalValues)]  ], 0.25*Quartile)

      # if evalValues is NA use a quartile (for vectors)
      if (length(evalNames[is.na(evalValues)])>1)
        evalData[evalNames[is.na(evalValues)]] <- apply(fit$x[,evalNames[is.na(evalValues)]  ], 2, quantile, 0.25*Quartile)


      # correct medians of interaction terms
      for (i in 1:length(iterms)){
        tmp <- items[items[,2]==iterms[i],1]
        if (length(tmp)>0)
          evalData[iterms[i]] <- prod(evalData[tmp])
      }

    return(evalData)
    }

# <-----------------------------------------------------------------------------
# (function starts here)



    RobVar <- fit$var
    if (is.null(RobVar)) RobVar <- vcov(fit)

    ShiftIntercept <- 0
    if (   length(grep("^Intercept$",  names(fit$coef)))>0 | length(grep("^\\(Intercept\\)$", names(fit$coef)))>0 ) ShiftIntercept <- 1

#    if (ShiftIntercept == 0) warning("No intercept found.")



    # list of interaction terms
    iterms <- grep('\\*', names(fit$coef))
    # list all variables which appear in an interaction
    items <- matrix(0, nrow=0, ncol=2)

    # list all variables which appear in an interaction
    for(i in 1:length(coef(fit))){
      itemsinteractions <- unique(c( grep(paste(names(fit$coef)[i], '\\*', sep=' '), names(fit$coef)), grep(paste('\\*', names(fit$coef)[i] , sep=' '), names(fit$coef)) ))
      if (length(itemsinteractions)>0){
        items <-  rbind(items,
                cbind(i, (itemsinteractions)) )}
    }

  # 1. Use median and medianspecial data to set values of evalData.
  #    Copy evalData to evalData1and evalData2. Overwrite all changed
  #    variables.

    # vector with medians (interaction-variables are corrected later)
    if (ShiftIntercept>0){
    evalData <- apply(cbind(Intercept=rep(1,dim(fit$x)[1]),
                            fit$x), 2, median)
    } else{
    evalData <- apply(fit$x, 2, median)

    }


    # replace names with =  AND SET VALUES TO 0 !!!!!

    fixedeffects <- grep("=", names(fit$coef))
    colnames(RobVar) <- gsub("=", "", colnames(RobVar))
    rownames(RobVar) <- gsub("=", "", rownames(RobVar))
    names(fit$coefficients) <- gsub("=", "", names(fit$coefficients))
    names(evalData) <- gsub("=", "", names(evalData))

    evalData[fixedeffects] <- 0



    #change evalData as given in medianspecial
    if (length(medianspecial)>1){
     evalData[medianspecial[seq(from=1, to=length(medianspecial), by=2)]] <- as.numeric(medianspecial[seq(from=2, to=length(medianspecial), by=2)])
    }




    gbefore   <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(1, 0), 1)

    gafter   <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(1, 1), 1)



    gbefore[which(names(gbefore) == treatment)] <- 1

  # replace interaction terms' names
        for (k in iterms){

            tmpl <- items[ items[,2]==k,1]
            tstr <- ''
            for (j in 1:length(tmpl)){
              tstr <- paste(tstr, tmpl[[j]], sep='x')
            }
            colnames(RobVar)[k] <- tstr
            rownames(RobVar)[k] <- tstr
            names(fit$coefficients)[k] <-tstr
            names(gbefore)[k] <- tstr
            names(gafter)[k] <- tstr
    }



    texE <- paste( LogitText(gafter), "-", LogitText(gbefore), sep="")

    tableDE <- deltaMethod(coef(fit), vcov.=RobVar, g=texE)
    rownames(tableDE) <- "Discrete Change"


    pvalue <- 2*(1-pnorm(abs(data.matrix(tableDE['Estimate']/tableDE['SE']))))
    colnames(pvalue)[1] <- 'pvalue'
    tableDE <- cbind(tableDE, pvalue)
class(tableDE) <- c(class(tableDE), "DiscreteEffects")
return(tableDE)
}






DiscreteEffectDiD <- function(fit, group, treatment, medianspecial=NA){

require(alr3)

    # check variable name: Intercept and "=" not allowed
    IsProperVariable <- function(i){
        if( names(fit$coef)[i]!='Intercept' && regexpr("=", names(fit$coef)[i])==-1 ){
        return(T)
        } else return(F)
      }


    # check whether variable with name vname has to be evaluated at
    # default value (=dval) or lower/upper exception value (=1, 2)

    evalexception <- function(vname, dval, lowup){
     if (length(evalspecial)>0  ){
      indx<- grep(vname, evalspecial[seq(from=1, to=length(evalspecial), by=3)])[1]
       if (length(indx)>0 & is.na(indx)==F){
           if (is.na( as.numeric(evalspecial[1+(indx-1)*3+lowup]))==F){
             return(as.numeric(evalspecial[1+(indx-1)*3+lowup]))
           }else{ stop("check evalspecial arguments")
                  return(dval)
                }
        }else{return(dval)}
     }else{return(dval)}
    }





    # construct Logit formula
    LogitText <- function(variables){
     txt <- ""
     for (i in 1:(length(variables)-1)){
       txt <- paste(txt, names(variables)[i], "*", variables[i], "+", sep=" ")
      }
     txt <- paste(txt, names(variables)[length(variables)], "*", variables[length(variables)],  sep=" ")

     txt <- paste( "exp(", txt, ")/ (1+ exp(",txt, "))", sep="")
     return(txt)
    }

    SetChangeVariables <- function(evalData, evalNames, evalValues, Quartile){

       evalValues <- as.numeric(evalValues)

      # values given by evalValues
      if (length(evalNames[is.na(evalValues)==F])>0)
        evalData[evalNames[is.na(evalValues)==F]] <- evalValues[is.na(evalValues)==F]

      # if evalValues is NA use a quartile (for one value)
      if (length(evalNames[is.na(evalValues)])==1)
        evalData[evalNames[is.na(evalValues)]] <- quantile(fit$x[,evalNames[is.na(evalValues)]  ], 0.25*Quartile)

      if (length(evalNames[is.na(evalValues)])>1)
      stop("Discrete Effects DiD Funktion �berpr�fen: Fehlende Werte werden auf Quartile = 1 gesetzt")

      # if evalValues is NA use a quartile (for vectors)
      if (length(evalNames[is.na(evalValues)])>1)
        evalData[evalNames[is.na(evalValues)]] <- apply(fit$x[,evalNames[is.na(evalValues)]  ], 2, quantile, 0.25*Quartile)


      # correct medians of interaction terms
      for (i in 1:length(iterms)){
        tmp <- items[items[,2]==iterms[i],1]
        if (length(tmp)>0)
          evalData[iterms[i]] <- prod(evalData[tmp])
      }

    return(evalData)
    }

# <-----------------------------------------------------------------------------
# (function starts here)



    RobVar <- fit$var
    if (is.null(RobVar)) RobVar <- vcov(fit)

    ShiftIntercept <- 0
    if (   length(grep("^Intercept$",  names(fit$coef)))>0 | length(grep("^\\(Intercept\\)$", names(fit$coef)))>0 ) ShiftIntercept <- 1

#    if (ShiftIntercept == 0) warning("No intercept found.")



    # list of interaction terms
    iterms <- grep('\\*', names(fit$coef))
    # list all variables which appear in an interaction
    items <- matrix(0, nrow=0, ncol=2)

    # list all variables which appear in an interaction
    for(i in 1:length(coef(fit))){
      itemsinteractions <- unique(c( grep(paste(names(fit$coef)[i], '\\*', sep=' '), names(fit$coef)), grep(paste('\\*', names(fit$coef)[i] , sep=' '), names(fit$coef)) ))
      if (length(itemsinteractions)>0){
        items <-  rbind(items,
                cbind(i, (itemsinteractions)) )}
    }

  # 1. Use median and medianspecial data to set values of evalData.
  #    Copy evalData to evalData1and evalData2. Overwrite all changed
  #    variables.

    # vector with medians (interaction-variables are corrected later)
    if (ShiftIntercept>0){
    evalData <- apply(cbind(Intercept=rep(1,dim(fit$x)[1]),
                            fit$x), 2, median)
    } else{
    evalData <- apply(fit$x, 2, median)

    }


    # replace names with =  AND SET VALUES TO 0 !!!!!

    fixedeffects <- grep("=", names(fit$coef))
    colnames(RobVar) <- gsub("=", "", colnames(RobVar))
    rownames(RobVar) <- gsub("=", "", rownames(RobVar))
    names(fit$coefficients) <- gsub("=", "", names(fit$coefficients))
    names(evalData) <- gsub("=", "", names(evalData))

    evalData[fixedeffects] <- 0



    #change evalData as given in medianspecial
    if (length(medianspecial)>1){
     evalData[medianspecial[seq(from=1, to=length(medianspecial), by=2)]] <- as.numeric(medianspecial[seq(from=2, to=length(medianspecial), by=2)])
    }




    gbefore   <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(1, 0), 1)

    gafter   <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(1, 1), 1)

    cbefore   <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(0, 0), 1)

    cafter    <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(0, 1), 1)


  # replace interaction terms' names
        for (k in iterms){

            tmpl <- items[ items[,2]==k,1]
            tstr <- ''
            for (j in 1:length(tmpl)){
              tstr <- paste(tstr, tmpl[[j]], sep='x')
            }
            colnames(RobVar)[k] <- tstr
            rownames(RobVar)[k] <- tstr
            names(fit$coefficients)[k] <-tstr
            names(gbefore)[k] <- tstr
            names(gafter)[k] <- tstr
            names(cbefore)[k] <- tstr
            names(cafter)[k] <- tstr
    }



    texE <- paste( "(", LogitText(gafter), "-", LogitText(gbefore), ") - (", LogitText(cafter), "-", LogitText(cbefore), ")", sep="")
    tableDE <- deltaMethod(coef(fit), vcov.=RobVar, g=texE)
    rownames(tableDE) <- "Discrete Change"


    pvalue <- 2*(1-pnorm(abs(data.matrix(tableDE['Estimate']/tableDE['SE']))))
    colnames(pvalue)[1] <- 'pvalue'
    tableDE <- cbind(tableDE, pvalue)
class(tableDE) <- c(class(tableDE), "DiscreteEffects")
return(tableDE)
}





#! Difference in Difference Select Bounds






DiscreteEffectDiDBounds <- function(fit, group, treatment, bounds.group, bounds.treatment, medianspecial=NA){

require(alr3)

    # check variable name: Intercept and "=" not allowed
    IsProperVariable <- function(i){
        if( names(fit$coef)[i]!='Intercept' && regexpr("=", names(fit$coef)[i])==-1 ){
        return(T)
        } else return(F)
      }


    # check whether variable with name vname has to be evaluated at
    # default value (=dval) or lower/upper exception value (=1, 2)

    evalexception <- function(vname, dval, lowup){
     if (length(evalspecial)>0  ){
      indx<- grep(vname, evalspecial[seq(from=1, to=length(evalspecial), by=3)])[1]
       if (length(indx)>0 & is.na(indx)==F){
           if (is.na( as.numeric(evalspecial[1+(indx-1)*3+lowup]))==F){
             return(as.numeric(evalspecial[1+(indx-1)*3+lowup]))
           }else{ stop("check evalspecial arguments")
                  return(dval)
                }
        }else{return(dval)}
     }else{return(dval)}
    }





    # construct Logit formula
    LogitText <- function(variables){
     txt <- ""
     for (i in 1:(length(variables)-1)){
       txt <- paste(txt, names(variables)[i], "*", variables[i], "+", sep=" ")
      }
     txt <- paste(txt, names(variables)[length(variables)], "*", variables[length(variables)],  sep=" ")

     txt <- paste( "exp(", txt, ")/ (1+ exp(",txt, "))", sep="")
     return(txt)
    }

    SetChangeVariables <- function(evalData, evalNames, evalValues, Quartile){

       evalValues <- as.numeric(evalValues)

      # values given by evalValues
      if (length(evalNames[is.na(evalValues)==F])>0)
        evalData[evalNames[is.na(evalValues)==F]] <- evalValues[is.na(evalValues)==F]

      # if evalValues is NA use a quartile (for one value)
      if (length(evalNames[is.na(evalValues)])==1)
        evalData[evalNames[is.na(evalValues)]] <- quantile(fit$x[,evalNames[is.na(evalValues)]  ], 0.25*Quartile)

      if (length(evalNames[is.na(evalValues)])>1)
      stop("Discrete Effects DiD Funktion �berpr�fen: Fehlende Werte werden auf Quartile = 1 gesetzt")

      # if evalValues is NA use a quartile (for vectors)
      if (length(evalNames[is.na(evalValues)])>1)
        evalData[evalNames[is.na(evalValues)]] <- apply(fit$x[,evalNames[is.na(evalValues)]  ], 2, quantile, 0.25*Quartile)


      # correct medians of interaction terms
      for (i in 1:length(iterms)){
        tmp <- items[items[,2]==iterms[i],1]
        if (length(tmp)>0)
          evalData[iterms[i]] <- prod(evalData[tmp])
      }

    return(evalData)
    }

# <-----------------------------------------------------------------------------
# (function starts here)



    RobVar <- fit$var
    if (is.null(RobVar)) RobVar <- vcov(fit)

    ShiftIntercept <- 0
    if (   length(grep("^Intercept$",  names(fit$coef)))>0 | length(grep("^\\(Intercept\\)$", names(fit$coef)))>0 ) ShiftIntercept <- 1

#    if (ShiftIntercept == 0) warning("No intercept found.")



    # list of interaction terms
    iterms <- grep('\\*', names(fit$coef))
    # list all variables which appear in an interaction
    items <- matrix(0, nrow=0, ncol=2)

    # list all variables which appear in an interaction
    for(i in 1:length(coef(fit))){
      itemsinteractions <- unique(c( grep(paste(names(fit$coef)[i], '\\*', sep=' '), names(fit$coef)), grep(paste('\\*', names(fit$coef)[i] , sep=' '), names(fit$coef)) ))
      if (length(itemsinteractions)>0){
        items <-  rbind(items,
                cbind(i, (itemsinteractions)) )}
    }

  # 1. Use median and medianspecial data to set values of evalData.
  #    Copy evalData to evalData1and evalData2. Overwrite all changed
  #    variables.

    # vector with medians (interaction-variables are corrected later)
    if (ShiftIntercept>0){
    evalData <- apply(cbind(Intercept=rep(1,dim(fit$x)[1]),
                            fit$x), 2, median)
    } else{
    evalData <- apply(fit$x, 2, median)

    }


    # replace names with =  AND SET VALUES TO 0 !!!!!

    fixedeffects <- grep("=", names(fit$coef))
    colnames(RobVar) <- gsub("=", "", colnames(RobVar))
    rownames(RobVar) <- gsub("=", "", rownames(RobVar))
    names(fit$coefficients) <- gsub("=", "", names(fit$coefficients))
    names(evalData) <- gsub("=", "", names(evalData))

    evalData[fixedeffects] <- 0



    #change evalData as given in medianspecial
    if (length(medianspecial)>1){
     evalData[medianspecial[seq(from=1, to=length(medianspecial), by=2)]] <- as.numeric(medianspecial[seq(from=2, to=length(medianspecial), by=2)])
    }




    gbefore   <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(bounds.group[2], bounds.treatment[1]), 1)

    gafter   <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(bounds.group[2], bounds.treatment[2]), 1)

    cbefore   <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(bounds.group[1], bounds.treatment[1]), 1)

    cafter    <- SetChangeVariables(evalData,
                                    c(group, treatment),
                                    c(bounds.group[1], bounds.treatment[2]), 1)


  # replace interaction terms' names
        for (k in iterms){

            tmpl <- items[ items[,2]==k,1]
            tstr <- ''
            for (j in 1:length(tmpl)){
              tstr <- paste(tstr, tmpl[[j]], sep='x')
            }
            colnames(RobVar)[k] <- tstr
            rownames(RobVar)[k] <- tstr
            names(fit$coefficients)[k] <-tstr
            names(gbefore)[k] <- tstr
            names(gafter)[k] <- tstr
            names(cbefore)[k] <- tstr
            names(cafter)[k] <- tstr
    }



    texE <- paste( "(", LogitText(gafter), "-", LogitText(gbefore), ") - (", LogitText(cafter), "-", LogitText(cbefore), ")", sep="")
    tableDE <- deltaMethod(coef(fit), vcov.=RobVar, g=texE)
    rownames(tableDE) <- "Discrete Change"


    pvalue <- 2*(1-pnorm(abs(data.matrix(tableDE['Estimate']/tableDE['SE']))))
    colnames(pvalue)[1] <- 'pvalue'
    tableDE <- cbind(tableDE, pvalue)
class(tableDE) <- c(class(tableDE), "DiscreteEffects")
return(tableDE)
}
